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We construct an accurate estimate for the root mean square force error of the particle-particle-parti- 
cle-mesh (P 3 M) algorithm by extending a single particle pair error measure which has been given by 
Hockney and Eastwood. We also derive an easy-to-use analytic approximation to the error formula. 
This allows a straightforward and precise determination of the optimal splitting parameter (as a 
function of system specifications and P 3 M parameters) and hence knowledge of the force accuracy 
prior to the actual simulation. The high quality of the estimate is demonstrated in several examples. 



o 



i 

C 

O 

o 



> 
o 
o 

o 

oo 
On 



T3 
C 

o 
o 



INTRODUCTION 

The combination of periodic boundary conditions and 
long range interactions is a frequent difficulty encoun- 
tered in computer simulations of physical systems, and 
the ingenious summation technique connected with the 
name of Ewalcu has become a standard instrument for 
tackling this problem. However, it has long been real- 
ized that for rather extensive simulations (involving, e.g., 
many particles) this approach is still too time consuming 
and various alternative methods have been invented. One 
particular class of these new algorithms owes its speedup 
to an inspired replacement of the Fourier transformation 
- which lies at the heart.oi the Ewald technique - by Fast 
Fourier (FFT) routinesacL 

In a recent publicationEl we presented a unified view 
of these FFT accelerated Ewald sums and carried out 
detailed accuracy measurements. However, since all these 
algorithms contain various free parameters, working at 
the maximally obtainable accuracy requires the user to 
tune them very carefully. This is straightforward if there 
exists a theoretical estimate of the errors involved - as 
is the case for the standard Ewald sumQ as well-, as for 
the so-called particle mesh Ewald (PME) methodQ - but 
rather tedious otherwise. 

In this article we present such an estimate for the 
root mean square error in the force of the so-called par- 
ticle-particle-particle-mesh (P 3 M) algorithm by extend- 
ing an erppr measure already derived by Hockney and 
Eastwooda and additionally provide an easy to use ana- 
lytical approximation to the somewhat unwieldy expres- 
sion comprising various sums. 



EWALD SUM AND P J M IN A NUTSHELL 

In this section we outline for reference purpose the 
most important formulas for the P 3 M method without 



much explanation, derivation or motivation. For the de= 
tails the reader is referred to the original P 3 M literatureo 
as well as our previous publicationlj. 

Consider a system of N particles with charges qi at 
positions in an overall neutral cubic simulation box of 
length L. Employing the Ewald sum, the force Fj on par- 
ticle i, which results from all interactions with the other 
charges (including all periodic images), can be written in 
the following way: 
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The so-called real space, Fourier space and dipole contri- 
butions are respectively given by 
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The prime on the second sum in Eqn. (^J) indicates 
that for i = j the term m = has to be omitted, 
erfc(r) := 2-7r~ 1 / 2 dt exp(— t 2 ) is the complementary 
error function and of course = — Tj. Further- 
more, the inverse length a is the splitting parameter of 
the Ewald sum, which controls the relative importance of 
the contributions coming from real and reciprocal space, 
the k- vectors are from the discrete set ^Z 3 and e' is the 
dielectric constant of the medium, which surrounds the 
cluster of simulation boxes as it tends (in a spherical way) 
towards an infinite system. In practice, the infinite sums 
in Eqns. (|2|, ||) are truncated by only taking into account 
distances which are smaller than some real space cutoff 
''max and wave vectors with a modulus smaller than some 
reciprocal space cutoff fc max . 
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The P 3 M method offers a fast way for an approximate 
computation of the reciprocal space contribution ([}]) . By 
mapping the system onto a mesh, the necessary Fourier 
transformations can be accomplished by Fast Fourier rou- 
tines. At the same time the simple Coulomb Green func- 
tion 4-7r/fc 2 is adjusted as to make the result of the mesh 
calculation most closely resemble the continuum solution. 

The first step, i.e., generating the mesh based charge 
density pm (defined at the mesh points r p ), is carried out 
with the help of a charge assignment function W: 
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Hereby the sum extends over the complete mesh M. 

Although the presented formulas d^-pj]) look some- 
what complicated, it is rather easy to implement them 
step by step. Furthermore, due to the replacement of 
the Fourier transforms by FFT routines (see Eqn. (0)), 
the algorithm is not only fast but its CPU time shows a 
favorable scaling with particle number: If the real space 
cutoff r max is chosen small enough (so that the real space 
contribution (]|) can be calculated in order A), the com- 
plete algorithm is essentially of order A log A. 



Here h is the mesh spacing, and the number of mesh 
points Am = L/h along each direction should preferably 
be a power of two, since in this case the FFT is most 
efficient. The charge assignment function is classified ac- 
cording to its order P, i.e. between how many grid points 
- per coordinate direction - each charge is distributed. In 
the P 3 M method introduced by Hockney and EastwoodB 
its Fourier transform is 
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In a second step the mesh based electric field E(r p ) 
is calculated. Basically, the electric field is the deriva- 
tive of the electrostatic potential, but there exist sev- 
eral alternatives for implementing the differentiation on 
a lattices. In this article we will restrict to the case of ik- 
differentiation, which works by multiplying the Fourier 
transformed potential with ik. In this case E(r p ) can be 
written as 



E(r„) = FFT 



-ik x FFT [ PM ] x G opt 



(r P ) (7) 



In words, E(r p ) is the backward finite Fourier transform of 
the product of — zk, the forward finite Fourier transform 
of the mesh based charge density pm and the so-called 
optimal influence function G Q pt, given by 
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Here D(k) is the Fourier transform of the employed dif- 
ferentiation operator, which is simply zk in our case. Fi- 
nally, one arrives at the force on particle i, i.e. the re- 
placement of Eqn. (||): 



SCALING OF THE RMS FORCE ERROR 

In this section we address the dependence of the root 
mean square error in the force on the number of charged 
particles and their valence. Since the assumptions and 
arguments involved are of a rather general nature, the 
result is not specific to a certain kind of Ewald method. 

We define the rms error in the force to be 



AF := 
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where Fj is the force on particle i calculated by the al- 
gorithm under investigation and Ff xa is the exact force 
on that particle. Note that this is by no means the only 
interesting measure of accuracy. However, it is the only 
one which is considered in this article. 

We now assume that the error in the force on particle 
i can be written as 



AFi = q % ^2 Qj Xi 



(13) 



The idea behind this ansatz is that - just as it is true for 
Fj - the error in Fj originates from the N — 1 interac- 
tions of particle i with the other charged particles, and 
each contribution should be proportional to the product 
of the two charges involved. The vector \ij gives the di- 
rection and magnitude of this error for two unit charges 
and depends on their separation and orientation as well 
as on the specific algorithm used for calculating the elec- 
trostatic forces. For this term we further assume 



(Xij ■ Xik) = S jk (xlj) = : Sjk X 



(14) 



where averaging over the particle configurations is de- 
noted by the angular brackets. The underlying assump- 
tion that contributions from different particles are un- 
corrected is certainly not always true (think, e.g., of 
highly ordered or strongly inhomogeneous particle dis- 
tributions), but it is sensible for random systems. Obvi- 
ously, the term (xfj) ~ the mean square force error for 
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two unit charges - can no longer depend on i and j and 
is thus written as x 2 - Using Eqns. (^Ju]), it follows 

((AF 4 ) 2 ) = q 2 ]T £ qMXij ■ X lk ) « X 2 Q 2 (15) 
where the important quantity Q 2 is defined as 
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Not all particles necessarily have the same charge. 
More specifically, let there be P subsets N p , defined by 



the condition that all 



particles from the subset N p 



have the same charge c v . If |N P | 3> 1, the law of large 
numbers and Eqn. (n5|) gives 



clx 2 Q 2 
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i.e., the arithmetic mean of the (AF^) 2 for all particles 
i £ N p can be approximated by the ensemble average for 
one particle from N p . In the case where all \N P \ are large, 
it follows 

»=1 p=l \ 1 Pl iGNp 
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Inserting this into Eqn. ( |12[ ) gives the final relation 
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Thus, the scaling of the rms error in the force with parti- 
cle number and valence is given by the factor Q 2 N^ 1 ^ 2 , 
whereas the prefactor x ~ which cannot be obtained 
by such simple arguments - contains the details of the 
method. Indeed, the estimates for the real and recipro- 
cal space error of the standard Ewald sumla as well as 
the estimate for the reciprocal space error of the PME 
methodQ are exactly of the form (|ll]). Note that any in- 
formation on the valence distribution enters only through 
the value of Q 2 . 



THE ERROR MEASURE OF HOCKNEY AND 
EASTWOOD 

The most interesting ingredient of the P 3 M method is 
the optimal influence function from Eqn. (|sj) . It is con- 
structed such that the result of the mesh calculation is 
as close as possible to the solution of the original con- 
tinuum problem. Of course, this is only realizable in a 
quantitative way, if the notion of "as close as possible" is 



stated more precisely. Hockney and Eastwood define the 
following measure of the error involved in a P 3 M calcu- 
lation: 



Q '- = hj h3 d3ri J L3 d3r t F(r; ri } _ R(r) ] : 
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F(r; ri) is the Fourier space contribution of the force be- 
tween two unit charges at positions ri and ri + r as cal- 
culated by the P 3 M method (note that due to broken 
rotational and translational symmetry this does in fact 
depend on the coordinates of both particles), and R(r) 
is the corresponding exact reference force (whose Fourier 
transform is just Eqn. (^)). The inner integral over r 
scans all particle separations, whereas the outer integral 
over ri averages over all possible locations of the first par- 
ticle within a mesh cell. Obviously, up to a factor L~ 3 
this expression is just the mean square error in the force 
for two unit charges, in other words, the quantity x 2 from 
Eqn. (|lj). This provides a link between the rms error of 
an N particle system and the error Q from Hockney and 
Eastwood: Using Eqn. (fl9|) one obtains 
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Technically spoken, Q is a functional of the influence 
function, and by setting the functional derivative of Q 
with respect to G to zero Hockney and Eastwood were 
able to derive the optimal influence function from Eqn. 
(||). However, it is most important to realize that they 
also provide a closed expression for the corresponding 
"optimal error" Q opt = Q[G op t\: 
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The outer sum extends over all k- vectors of the Fourier 
transformed mesh M, and the star denotes complex 
conjugation. Once again, in the special case of ik- 
differentiation one has D(k) = ik. 

Admittedly, Eqn. (p3) looks rather complicated. Still, 
in combination with Eqn. ( ^l|) it gives the rms force error 
of the P 3 M method (or - more precisely - of its Fourier 
space contribution)! (After all, the computation of Q pt 
and that of G op t are quite similar.) We would like to 
emphasize that the formula ( p2| ) for the optimal Q- value 
(just like the one for the optimal influence function (^)) 
is of a very general nature: It does also work for dif- 
ferent charge assignment functions, reference forces or 
any differentiation scheme which can be expressed by an 
operator D(k), in particular for all the finitC|-flifference 
schemes presented in our previous publicationO. 
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The correspondingrms error in the force from the real 
space contribution (g) has been derived by Kolafa and 
PerramQ and we want to provide it here for reference 
purpose: 



with the function defined as 
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In a second step the aliasing sum in the denominator 
of Eqn. (|2^) is evaluated analytically by exploiting the 
following partial fraction expansionB: 



sin 2 {x) = [x + rrm)' 



x e M\ttZ (27) 



With these two estimates at hand it is easy to deter- 
mine the optimal value of the splitting parameter a via a 
stand-alone program, which takes the relevant system pa- 
rameters (N, Q 2 , L) and specifications of the algorithm 

(r max , A^m, P) as its input. If real and reciprocal space Differentiating this expression 2P - 2 times gives 
contribution to the error, AF l - r > and AF 1 - > respectively, 
are assumed to be statistically independent, the total er- 
ror is given by 
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(24) 



This quantity has to be minimized with respect to a. 
However, in most cases it is accurate enough to use the 
following approximation: Determine the value of a at 
which the real and reciprocal space contribution to the 
rms force error are equal. 



ANALYTIC APPROXIMATION 

Although the closed expression for the error from the 
last section is not really complicated, it is somewhat un- 
wieldy. A possible calculation of the optimal value of a 
using e.g. a bisection method needs several computations 
of Qopt and for each it is necessary to compute the in- 
ner aliasing sums and the outer sum over the k- vectors. 
Especially for large Fourier meshes this can be rather 
time consuming. Therefore we now derive an analytic 
approximation to this error estimate, which is essentially 
an expansion for small ha. We will restrict to the case of 
a cubic system and the same number Nm of mesh points 
along each direction. Also, we will only treat the. case 
of ik-differentiation (see Eqn. ([?])), since we foundcl that 
this is the most accurate method. However, our line of 
reasoning can be extended to more general cases. 

We start our treatment of Eqn. (g2|) by observing that 
two of the three sums over Z 3 contain R with its expo- 
nential factor exp(— fc 2 /4a 2 ). Since near the boundary 
of M its value is roughly given by exp(— (ir/ha) 2 ), R is 
strongly damped outside M if ha is small. Thus, it is a 
good approximation to retain only the term with m = 
in these two sums. Inserting D(k) = ik, the Fourier 
transform of the charge assignment function from Eqn. 
(||) and using the fact that sin 2 (x + mr) — sin 2 (x) for 
integral n, one obtains 
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This equation leads to a closed expression for the function 
/( p ) from Eqn. <M): 
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Unfortunately the sum over M is still too complicated to 
perform, so some further approximations are necessary. 
We choose the following way: is expanded in a Tay- 
lor series up to order AP — 2. Since (i) f( p ' is an even 
function, (ii) f( p \0) = 1 and (Hi) the lowest nontrivial 
term is of order x 2P , this expansion can be written as 



f^(x)^fP(x) -.= i-x- ^ *■:.■>■ 
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The coefficients c m are easily determined with the help 
of any mathematical computer program capable of sym- 
bolic algebra and are listed in Table |. The term in curly 
brackets from Eqn. (25) can now be approximated like 
this: 

i _ f{p)(h^)f{p)(hft)f(p)(h]l) Rj 



(31) 

The product of the three functions is computed by 
multiplying their Taylor expansions term by term, but 
the results are only retained up to the truncation order 
4P — 2. Note that the first neglected cross term would 
be of order 4P. 

For symmetry reasons it is clear that all three terms 
in the last line of (|3l]) contribute in the same way to 
the value of the sum in Eqn. p5|), therefore it suffices to 
choose one of them - e.g., the z-term - and multiply the 

(p) 

result by 3. Together with the definition of fj, from 
Eqn. (pG) this leads to 
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Finally, the sum is replaced by an integral via 
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If one extends the range of integration to M 3 and changes 
to spherical polar coordinates, the remaining angular and 
radial integrals can be performed with the help of 
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where (2n — 1)!! = 1 ■ 3 ■ 5 • ■ ■ (2n — 1). Collecting all parts 
together gives 
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with the abbreviation 
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Combining this with Eqn. ( ]21| ) results in the following 
analytical approximation to the Fourier space contribu- 
tion to the rms force error of the P 3 M algorithm: 
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The exact expansion coefficients a m (which are rational 
numbers) are listed in Tab le p . 

Let us repeat that Eqn. (|38|) was derived under the ex- 
plicit assumption that ha is small. Both, the restriction 
to the term m = for two sums in Eqn. ( |22| ) as well 
as the expansion of the function /( p ) from Eqns. ( |26| , p9|) 
in powers of ha can become questionable if ha becomes 
large. However, in this case it is still safe to go back to 
the original error estimate, i.e. the combination of Eqns. 
(HI) and (||). 



NUMERICAL TEST 

In this section we demonstrate the accuracy of the P 3 M 
error estimates by comparing their predictions with the 
exact rms force error AF from Eqn. (||) - calculated 
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FIG. 1. The rms error AF (solid lines) for the system 
of 100 randomly distributed charges is calculated for the 
ik-differentiated P 3 M method with Nyi = 32 mesh points 
and real space cutoff r max = 4. From top to bottom the or- 
der of the charge assignment function is increased from 1 to 7. 
The dotted lines are the corresponding full estimates (using 
Eqns. (plp2|)) for the Fourier space contribution to AF. 
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FIG. 2. Same plot as in Fig. [l], but here the dotted lines 
are the analytical estimates from Eqn. Note that for 

large P the error AF is overestimated at large values of a. 



for a specific random system. Hereby the exact forces Fj 
needed for computing AF are obtained by a well con- 
verged standard Ewald sum, and for the test system we 
choose the |One described in Appendix D of our previous 
publication^: 100 particles randomly distributed within 
a cubic box of length L = 10, half of them carry a pos- 
itive, the other half a negative unit charge. Our unit 
conventions are as followaj: lengths are measured in C 
and charges in C. Hence the unit of force is C 2 /£ 2 . We 
will refer to the estimate which emerges from combining 
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FIG. 3. AF (solid lines) for the same system as in the 
previous figures is calculated for the ik-differentiated P 3 M 
method with charge assignment order P = 3 and real space 
cutoff r max = 4. From top to bottom the number of mesh 
points varies like 4, 8, 16, 32, 64, 128. The dotted lines are 
the corresponding full error estimates. 



FIG. 5. Test of the Q 2 N~ 1/2 scaling of AF. The three 
solid lines show the rms force error AF for systems, which 
differ in their (Q , N) values. From top to bottom they are 
characterized by (10000,400), (1000,200) and (100, 100) (the 
last system is the same as the one in Figs, hl-kj) . The dotted 
curves are the corresponding full error estimates. 




FIG. 4. Same plot as in Fig. but here the dotted lines 
are the analytical estimates from Eqn. ([^8|). At small N^, 
which corresponds to large values of h = L/Nm, the error 
formula overestimates AF at large values of a. 

Eqns. ( [2l|) and (^2|) as the full estimate and to Eqn. (pq) 
as the analytical approximation. 

In a first example we fix the number of mesh points 



to N 



M 



32 and the real space cutoff to 



The charge assignment order varies from P = 1 through 
P = 7. In Fig. |l| the resulting curves for the rms force 
error AF are plotted together with the full error estimate 
and in Fig. ||the same is done for the analytical approxi- 
mation. It can be seen very clearly that the full estimate 
accurately predicts the Fourier space contribution to AF 



for all values of a -and P. Since the real space contribu- 
tion is also knownD (see also Eqn. (|23|)), this permits an 
easy determination of the optimal value of the splitting 
parameter a. The analytical approximation is almost as 
accurate as the full formula, however, for large P it di- 
minishes in accuracy if a gets large. This is due to the 
fact that Eqn. (3^) was derived under the assumption 
that ha is small. Note that the expansion coefficients 
ain ^ needed in Eqn. (|3^) strongly increase with increas- 
ing m if P gets larger (see Table ||) . Still, both estimates 
are useful for determining the optimal operation point, 
and Eqn. ( J38| ) can be calculated much faster than the 
sums from Eqn. (p2|). 

Next we study at fixed charge assignment order P = 3 
different mesh sizes h — L/Nm by investigating iV*M S 
{4, 8, 16, 32, 64, 128}. Figs. | and | show AF in compar- 
ison with the full and the approximated error estimate 
respectively. Again, it can be seen that the former very 
accurately gives the Fourier space contribution to AF. 
As expected, the analytical approximation has problems 
at small Nm (since this results in large h) , but neverthe- 
less it is very useful otherwise. Essentially, one has to 
check the value of ha: If this is of the order unity or even 
larger, care is called for. Note that in Fig. [| the value of 
ha is approximately 1 at the points where the analytical 
approximation starts to deviate from the true curve. 

In a next step we want to demonstrate that the scaling 
of the rms force error with particle number and valence 
distribution is in fact correctly given by AF cx Q 2 N^ X / 2 . 
To this end we investigate three systems which differ only 
in the values of Q 2 and N. The first system is the same 
as the one investigated so far in Figs. A second 

system contains 200 particles, namely, 50 monovalent and 
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FIG. 6. Measured relative frequency p(r/L) for the scaled 
minimum image separation r/L between two charges of the 
polyelectrolyte system described in the text (solid line). For 
comparison, the probability density of r/L for a random sys- 
tem is also shown (dashed line). 



50 trivalent pairs. Finally, a third system contains 400 
particles: 50 pairs with charge ±1, 100 pairs with charge 
±5 and 50 pairs with charge ±7. Hence, their (Q 2 ,N) 
values are respectively given by (100,100), (1000,200) 
and (10000,400), and the ratio of their scaling prefactors 
is thus 1 : -s/50 : 50. In Fig. || this is clearly visible in 
the constant shift of the three curves with respect to one 
another (note that the vertical scale is logarithmic). Also 
shown is the full error estimate, which again predicts the 
Fourier space contribution to AF in all three cases very 
precisely. 



APPLICATION TO AN INHOMOGENEOUS 
POLYELECTROLYTE SYSTEM 

So far we have only used homogeneous random sys- 
tems for testing the error estimates. However, this does 
not necessarily reflect the situation encountered in all 
computer experiments. In the present section we want 
to show that deviations from a random distribution, as 
they frequently occur in charged systems, in fact have a 
noticeable influence on the rms error. 

We will use a typical system from our own research 
to demonstrate rtbis effect: a simple model of a polyelec- 
trolyte solutionBlilJ. 106 Lennard- Jones particles were 
joined (by some bonding potential) to build up a poly- 
mer chain. Every third "monomer" was monovalcntly 
charged and 8 such chains together with 96 trivalent and 
oppositely charged counterions, which make the complete 
system electrically neutral, were put in a cubic simulation 
box of length L ss 179. The system was brought into the 
canonical state by means of a Molecular Dynamics sim- 
ulation and a Langevin thermostat. 



FIG. 7. rms force error AF for the polyelectrolyte system 
described in the text (solid line). The dotted curve is the 
full P 3 M estimate for the reciprocal space contribution, the 
dashed curve is the estimate for the real space contribution. 
Note that due to the strong inhomogeneities in the charge 
distribution (see Fig. ^) both estimates show systematic de- 
viations from the true error curve. 



Under certain circumstances (e.g., at sufficiently low 
temperature and the appropriate density range) such 
polyelectrolyte chains collapse, and this happened to 
the described system. The chain sizes shrunk to much 
smaller values than for comparable neutral polymers and 
90% of the counterions were condensed within a distance 
of only two Lennard- Jones radii from the nearest chain. 

The various phenomena leading to this transition, the 
influence of the system parameters or the dynamics are 
only a few of the interesting physical questions. How- 
ever, the only thing which concerns us here is the fact 
that after the transition the system has developed local 
inhomogeneities. To demonstrate this we have plotted 
in Fig. U the measured relative frequency p(r/L) of the 
scaled minimum image distance r/L S [0; \/3/2] between 
two charges (we did not distinguish between different va- 
lences) and compared this to the probability density of 
r/L for a random system. The differences are indeed very 
pronounced. Apart from the more complicated structure 
of the measured curve, note in particular that small sep- 
arations are more frequent at the expense of larger ones. 

For this system we calculated the rms force error AF 
and the corresponding estimates for the real and recip- 
rocal space contribution (Eqns. (p3|;pl 22)). Since the 



simulation box comprises 288 monovalent and 96 triva- 
lent charges, we have N = 384 and Q = 1152. The 
result is shown in Fig. 0. No longer do the estimates 
correctly predict the two branches of AF. Rather, at 
small values of a the algorithm gives better results than 
expected from Eqn. (|23|), while at large values of a the 
estimated AF is smaller than the actual one. However, 
this trend can be explained qualitatively in the following 
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way: At small values of a the force (and also its error) 
is dominated by the real space contribution from Eqn. 
(||). This error originates from neglected contributions 
beyond the real space cutoff r max . In the case of Fig. [?] 
we used r max /L 0.243 and from Fig. || it can be seen 
that there are more particles within this cutoff (and thus 
less beyond) than in the case of a random particle dis- 
tribution, which should lead to an enhanced real space 
accuracy. On the other hand, we demonstrated previ- 
ously (see Fig. 7 of our recent publications) that the rms 
error of the P 3 M method strongly increases with decreas- 
ing minimum image distance. The general shift towards 
smaller r, which can be observed in the polyelectrolyte 
system, should thus lead to an enlarged reciprocal space 
error. Observable effects will occur at large values of a, 
where this contribution to AF dominates. 

Although the systematic deviations of the error esti- 
mates from the true curve are easily detectable, they 
are less dramatic than one could have guessed from a 
brief look at Fig. 0. The optimal splitting parameter 
from Fig. is given by a op t ~ 0.0715 with a correspond- 
ing AF op t ~ 1.2 x 10~ 6 , while the intersection point of 
real and reciprocal space estimate occurs at a ~ 0.0740, 
which predicts an error of AF ss 9.3 x 10~ 7 . If the esti- 
mated value of a had been used, this would result in an 
error of AF sa 1.5 x 10~ 6 , which is roughly 25% larger 
than at the optimal value of a. If such a safety margin is 
considered already at the beginning, the a priori deter- 
mination of the "optimal" value of a by means of Eqns. 
(|T]j2^j2^) together with an a posteriori validity check is 
still a good approach. 

In any case, if one knows or at least has reasons to 
suspect that the investigated system is susceptible to the 
development of inhomogeneities, one should always be 
aware of a potential failure of the presented error formu- 
las. In case of doubt, some simple numerical tests - like, 
e.g., the ones which we have performed here - are neither 
out of place nor costly. 

CONCLUSIONS 

We have presented an accurate estimate of the root 
mean square error in the force involved in a P 3 M cal- 
culation and additionally derived an easy-to-use analytic 
approximation for the special case of zk-diffcrentiation. 
Together with the existing estimate for the real space 
contribution to the rms force error this permits a deter- 
mination of the optimal tuning parameter a and provides 
information on the accuracy which is to be expected. It 
thus guarantees that one does not waste accuracy which 
could otherwise be achieved at the same computational 
effort. Stated the other way around: It prevents one 
from spending more computational effort for a desired 
algorithmic accuracy than actually necessary. 

In our previous publication^] we showed that the ik- 
differentiated P 3 M method is the most accurate algo- 



rithm for the FFT accelerated Ewald sum. In combina- 
tion with the error estimates this gives a "package" for 
calculating electrostatic interactions in periodic bound- 
ary conditions which is easy to use, very precise, produces 
known and controllable errors and can thus be optimally 
tuned in advance. If o«e wishes to rely on the discrete dif- 
ferentiation operatorsO, the full error formula ( ^l| , p2| ) will 
still work. It is. only in the case of the analytic differen- 
tiation schemeQ that none of the present error estimates 
is applicable. 

As a last point, we stressed several times that the va- 
lidity of the error estimates is subject to some additional 
requirements, concerning e.g. the homogeneity of the sys- 
tem. Nevertheless, it should be obvious that in any case 
a consultation of the error formulas - perhaps only for a 
first starting point - is superior to guessing the parame- 
ters or the use of a- values which were historically handed 
down. The difficulties of the analytical approximation 
( j38| ) at large values of ha are not a serious problem any- 
way, since in case of doubt one can always go back to the 
full estimate. 
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